% Copyright (C) 2009,2010,2011,2012  Marco Restelli
%
% This file is part of:
%   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
%
% LDGH is free software: you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% LDGH is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
% License for more details.
%
% You should have received a copy of the GNU General Public License
% along with LDGH. If not, see <http://www.gnu.org/licenses/>.
%
% author: Marco Restelli                   <marco.restelli@gmail.com>


function PV_dgcomp_series(fname,itime,N,k,l,outfile,varargin)
% PV_dgcomp_series(fname,itime,N,k,l,outfile)
% PV_dgcomp_series(fname,itime,N,k,l,outfile,select)
%
% Call PV_dgcomp for a series of files and prepare a .pvd file to
% collect the generated vtk files.
%
% fname nust contain a '%' to indicate the subdomain number and a '$'
% to indicate the suffix 'grid', 'base' or 'res-XXXX'.
%
% The optional argument 'select' can be used as follows:
% 0 : no postprocessing
% 1 : total values, i.e. deviations + IO reference state
% 2 : total values; include also all the additional diagnostics
% 3 : primitive variables deviations from IO reference

 % input arguments
 total = 0;
 alldiags = 0;
 primdevs = 0;
 if(nargin>6)
   switch varargin{1}
    case 1
     total = 1;
    case 2
     total = 1;
     alldiags = 1;
    case 3
     primdevs = 1;
    otherwise
     error("Invalid value for 'select'");
   endswitch
 end

 % Summary file
 fnameNX = [outfile,'.pvd'];
 fid = fopen(fnameNX,'w');
 
 line = '<?xml version="1.0"?>';
 fprintf(fid,'%s\n',line);
 line = '<VTKFile type="Collection" version="0.1" byte_order="BigEndian">';
 fprintf(fid,'%s\n',line);
 line = ' <Collection>';
 fprintf(fid,'%s\n',line);

 % Part loop
 Npos = find(fname=='%');
 for j=1:length(itime)
   for i=0:N
     fnameN = [fname(1:Npos-1),num2str(i,'%.3i'),fname(Npos+1:end)];
     Xpos = find(fnameN=='$');
     fnameNX = [fnameN(1:Xpos-1),'grid',fnameN(Xpos+1:end)];
     grid = load(fnameNX);
     fnameNX = [fnameN(1:Xpos-1),'base',fnameN(Xpos+1:end)];
     base = load(fnameNX);
     fnameNX = [fnameN(1:Xpos-1),'res-',num2str(itime(j),'%.4i'), ...
                fnameN(Xpos+1:end)];
     res = load(fnameNX);

     if(total==1)
       fnameNX = [fnameN(1:Xpos-1),'init',fnameN(Xpos+1:end)];
       ref = load(fnameNX);
       uuu = res.uuu + ref.uuu_ioref;
       uuu(1,:,:) = uuu(1,:,:) + permute(ref.atm_ref.rho,[3 1 2]);
       uuu(2,:,:) = uuu(2,:,:) + permute(ref.atm_ref.e,[3 1 2]);
     elseif(primdevs==1)
       fnameNX = [fnameN(1:Xpos-1),'init',fnameN(Xpos+1:end)];
       ref = load(fnameNX);
       [tmp,uuu,tmp,pitheta] = dgcomp_postp(res.uuu , base.base , ...
                              ref.uuu_ioref , ref.atm_ref , ref.phc);
       % This is a trick to compute pi and theta for the reference
       % state, so that we can plot the deviations.
       [tmp,tmp,tmp,io_pith] = dgcomp_postp(0*res.uuu , base.base , ...
                              ref.uuu_ioref , ref.atm_ref , ref.phc);
     else
       uuu = res.uuu;
     end

     fnameNX = [outfile,'-',num2str(itime(j)), ...
                        '-',num2str(i,'%.3i'),'.vtu'];
     if(alldiags)
       PV_dgcomp(grid.grid,base.base,uuu,k,l,fnameNX,res.test_name,res.add_diags);
     elseif(primdevs)
       PV_dgcomp(grid.grid,base.base,uuu,k,l,fnameNX,res.test_name,pitheta-io_pith);
     else
       PV_dgcomp(grid.grid,base.base,uuu,k,l,fnameNX,res.test_name);
     end
     line = ['  <DataSet timestep="', num2str(res.time), ...
             '" group="" part="', num2str(i,'%i'), ...
             '" file="',fnameNX,'"/>'];
     fprintf(fid,'%s\n',line);
   end
 end

 line = ' </Collection>';
 fprintf(fid,'%s\n',line);
 line = '</VTKFile>';
 fprintf(fid,'%s\n',line);

 fclose(fid);

return

